Multi-particle collision dynamics modeling of viscoelastic fluids 
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In order to investigate the rheological properties of viscoelastic fluids by mesoscopic hydrodynam- 
ics methods, we develop a multi-particle collision dynamics (MPC) model for a fluid of harmonic 
dumbbells. The algorithm consists of alternating streaming and collision steps. The advantage of 
the harmonic interactions is that the integration of the equations of motion in the streaming step 
can be performed analytically. Therefore, the algorithm is computationally as efficient as the orig- 
inal MPC algorithm for Newtonian fluids. The collision step is the same as in the original MPC 
method. All particles are confined between two solid walls moving oppositely, so that both steady 
and oscillatory shear flows can be investigated. Attractive wall potentials are applied to obtain a 
nearly uniform density everywhere in the simulation box. We find that both in steady and oscilla- 
tory shear flow, a boundary layer develops near the wall, with a higher velocity gradient than in the 
bulk. The thickness of this layer is proportional to the average dumbbell size. We determine the 
zero-shear viscosities as a function of the spring constant of the dumbbells and the mean free path. 
For very high shear rates, a very weak "shear thickening" behavior is observed. Moreover, storage 
and loss moduli are calculated in oscillatory shear, which show that the viscoelastic properties at 
low and moderate frequencies are consistent with a Maxwell fluid behavior. We compare our results 
with a kinetic theory of dumbbells in solution, and generally find good agreement. 

PACS numbers: 47.ll.-j, 83.60.Bc, 66.20.+d 



I. INTRODUCTION 

It is the characteristic feature of soft matter sys- 
tems that a macromolecular component of nano- to mi- 
crometer size is dispersed in a solvent of much smaller 
molecules. The mesoscopic length scale of the dispersed 
component implies that crystalline phases have a very 
small shear modulus - which roughly scales like the in- 
verse of the third power of the structural length scale 
- and that both crystalline and fluid phases are charac- 
terized by long structural relaxation times. Soft matter 
systems have therefore interesting dynamical properties, 
because the time scale of an external perturbation can 
easily become comparable with the intrinsic relaxation 
time of the dispersed macromolecules. 

One of the unique properties of soft matter is its vis- 
coelastic behavior [lj. Due to the long structural relax- 
ation time, the internal degrees of freedom cannot relax 
sufficiently fast in an oscillatory shear flow, so that there 
is some elastic restoring force which pushes the system 
back to its previous state. A very well studied example 
of viscoelastic fluids are polymer solutions and polymer 
melts Q, Si. 

In the case of polymer melts, the charac- 
teristic time scale is given by the reptation time, i.e. by 
the time it takes a chain to slide by its contour length 
along the tube formed by other polymer chains [3|. 
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In order to bridge the length- and time-scale gap be- 
tween the solvent and macromolecular or colloidal scales, 
several mesoscopic simulation techniques - such as the 
lattice-Boltzmann method, dissipative-particle dynamics 
(DPD), and multi-particle collision dynamics (MPC) - 
have been suggested in recent years, and are in the pro- 
cess of being developed further. The idea of all these 
methods is to strongly simplify the microscopic dynamics 
in order to gain computational efficiency, but at the same 
time to exactly satisfy the conservation laws of mass, 
momentum and energy, so that hydrodynamic behavior 
emerges naturally on larger length scales. 

We will focus here on the multi-particle collision dy- 
namics (MPC) techniquel, i, @, also called stochas- 
tic rotation dynamics 7] (SRD), originally developed for 
Newtonian fluids. This particle-based hydrodynamics 
method consists of alternating streaming and collision 
steps. In the streaming step, point particles move bal- 
listically. In the collision step, particles are sorted into 
the cells of a simple cubic (or square) lattice. All 
particles in a cell collide by a rotation of their ve- 
locities relative to the center-of-mass velocity around 
a random axis [J]. A random shift of the cell lattice 
is performed before each collision step in order to re- 
store Galilean invarianceQ- This method has been ap- 
plied very successfully to study the hydrodynamic be- 
havior of many complex fluids, such as polymer solu- 
tions in equilibrium [H, Q and flowfiol. [ill lip - colloidal 
dispersions] 13, 14, vesicle suspensions [15|, |l6(, and reac- 
tive fluids (17I. fl8j| . 

The viscoelastic behavior of polymer solutions leads 
to many unusual flow phenomena, such as shear- 
induced phase separation [H, [2(| HH, viscoelastic phase 
separation 22], and elastic turbulence (23|. A coarse- 
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grained description of viscoelastic fluids is necessary in 
order to obtain a detailed understanding of the role of 
elastic forces in such flow instabilities. 

However, there is a second level of complexity in soft 
matter system, in which a colloidal component is dis- 
persed in a solvent, which is itself a complex fluid. 
Examples are spherical or rod-like colloids dispersed 
in polymer solutions or melts, which are exposed to 
a shear flow^I HE HE H3, S|- Shear flow can in- 
duce particle aggregation and alignment in these sys- 
tems. This is important, for example, in the processing 
of nanocomposites (28j . 

The aim of this paper is therefore the development of 
a MPC algorithm, which is able to describe viscoelastic 
phenomena, but at the same time retains the computa- 
tional simplicity of standard MPC for Newtonian fluids, 
and thereby allows to take advantage of this mesoscale 
simulation for the investigation of flow instabilities as well 
as suspensions with viscoelastic solvents. We show that 
this goal can be achieved by replacing the point particles 
of standard MPC by harmonic dumbbells. In order to 
obtain a strong elastic contribution, we consider a fluid, 
which consists of dumbbells only. However, it is of course 
straightforward to mix dumbbells with a point-particle 
solvent. A similar idea has been suggested recently for 
DPD fluids (29|. 



II. THE MODEL 

A. Algorithm 

In our MPC model, we consider N p point particles of 
mass m, which are pairwise connected by a harmonic 
potential V(ri,r2) = ^K(ri — r 2 ) 2 to form dumbbells, 
where K is the spring constant. The center-of-mass 
position r? and velocity v? for each dumbbell i, with 
i = 1, 2, Np/2, are represented by 

= ^(r,:i + r i2 ) ; = -(v a + v l2 ) . (1) 

Here r$i, and Vji, denote the position and veloc- 
ity of the two point particles composing a dumbbell i, 
respectively. 

The MPC algorithm consists of two steps, streaming 
and collisions [il [E l30j| . In the streaming step, within a 
time interval h, the motion of all dumbbells is governed 
by Newton's equations of motion, 

dt 4 ' dt 1 ' 1 ' 

where m c — 2m is the mass of a dumbbell, and if is 
the total external force on dumbbell i. We consider only 
constant force fields. The center-of-mass positions and 
velocities of dumbbells are then given by a simple ballis- 
tic motion. The evolution of the relative coordinates of 



each dumbbell are determined by the harmonic interac- 
tion potential, so that 

r n (t + h) — r i2 (t + h) = Ai(t) cos(tu h) 

+B l (t) sm{uj Q h) ; (3) 
Vji(i + h) — Vi2(t + h) = -LUoA^t) sin(w /i) 

+u> Bi(t) cos(u h) , (4) 

with angular frequency luq = -\/2K/m. The vectors 
and Hi arc different for each time step, and are calculated 
from the relative positions and velocities of the point par- 
ticles of dumbbell i before the streaming step, 

Ai(t) = r«(f) - r l2 (t) ; Bi(t) = — (va(t) - v a (t)) . 

UJ 

(5) 

In the MPC algorithm described here, r c , v c , A and B 
are continuous variables, evolving in discrete increments 
of time. In the absence of she ar flow , the average l ength 
of the dumbbell is r (d) = V / P)^ 1 = y/ d k sT/K for a 
rf-dimensional system. 

In the collision step, the point particles are sorted into 
the cells of a cubic lattice with lattice constant a,Q. Multi- 
particle collisions are performed for all particles in a cell 
J, by the same SRD algorithm^ as for point particle flu- 
ids. The velocity of each particle relative to the center- 
of-mass velocity v cm , j of the cell is rotated around a ran- 
domly chosen axis by a fixed angle a, 

v'j (t + h)= v cm , j + 11(a) [ V j (t + h)- v cm , j] , (6) 
where IZ(a) is a stochastic rotation matrix, and 

Nj 

with Nj the number of particles within cell J. This step 
guarantees that each particle changes the direction as 
well as the magnitude of its velocity during the multi- 
particle collisions, while the local momentum and the ki- 
netic energy are conserved. Random shifts are applied in 
each direction, so that the Galilean invariance is ensured 
even in case of small mean free path[7l. l3l(. 

In order to describe Couette or oscillatory shear flow, 
the system is confined within two parallel hard walls in 
the y direction, which are moving oppositely along the x 
direction. Here, L Xl L y and L z are used to denote the 
dimension of the simulation box along the correspond- 
ing directions. For a steady shear flow, the shear rate is 
given by 7 = 2v Vfa n >x /L y , with Uwaii^ the x component of 
the velocity of the wall moving along the positive direc- 
tion. Periodic boundary conditions are applied in the x 
and z direction, bounce-back boundary condition in the 
y direction. The system is therefore divided into L x /ao 
and L z /ao cells in the x and z directions (parallel to the 
walls) , but L y / ao + 1 cells in the y direction because of the 
random shifts. At the walls, for collision cells which are 
not completely filled by particles, extra virtual point par- 
ticles are added to conserve the monomer number density, 
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p, defined by the average number of monomers per cellQ. 
In principle, the velocities of the virtual particles can be 
drawn from a Maxwell-Boltzmann distribution of average 
velocity equal to the wall velocity and variance \Jk B T /m, 
where k B T the bulk temperature. In the simulation code, 
it is not necessary to sample the velocity of virtual wall 
particles individually. A random vector from Maxwell- 
Boltzmann distribution with wall velocity and variance 
\J{p — n)k B T/uY is then used instead of the contribution 
of the entire virtual particles in the cell, where n is the 
number of real particles in that cell. For point particles, 
the combination of bounce-back boundary condition and 
virtual wall particles has been shown to guarantee no-slip 
boundary condition to a very good approximation^. 



B. Thermostats 

In order to keep the system temperature constant, var- 
ious thermostats can be employed. In the first case, 
the MPC method with collisions by stochastic rotations 
(MPC-SRD) of relative velocities is augmented by ve- 
locity rescaling. The simulation box is subdivided into 
Ly/ao layers parallel to the walls. In each layer, the new 
velocity of each particle j in cell J is obtained by 
rescaling the velocity relative to the center-of-mass ve- 
locity of that cell, 



Vj = V cm ,J + (Vj - V cm ,j) 



knT> 



(8) 



Here k B T' is calculated from the actual velocity distri- 
bution 



N. 
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JtElaycr j — 1 j£ layer 

(9) 

where Nj denotes the number of particles in cell J and 
-Wiayer the number of cells which contains particles within 
a layer. 

In the second case, the Anderson's thermostat ver- 
sion of MPC, denoted MPC- AT, is applied 0,11. This 
thermostat employs a different collision rule instead 
of Eq. ©. In the MPC-AT-a version of the algo- 
rithm (without angular momentum conservation, com- 
pare Sec. Ill Cl below). the new velocities of point particles 
in the collision step are assigned as[32| 



v cm , J 



N K r 



N, 



fe=i 



A' 



(10) 



Here vj an is a velocity chosen from the Maxwell- 
Boltzmann distribution and Nk the number of particles 
within cell K. Instead of energy conservation in MPC, 
the temperature is kept constant in MPC- AT. 



C. Angular Momentum Conservation 

The standard MPC algorithm as well as the Ander- 
son thermostat version do not conserve angular momen- 
tum. It has been shown recently[34j that this lack of 
angular-momentum conservation may lead to quantita- 
tive or even qualitative incorrect results, like non-physical 
torques in circular Couette flows. We therefore also con- 
sider the ang ular-momentum conserving modification of 
MPC-ATflliil, denoted MPC-AT+a. Here, the veloc- 
ities in the collision step are calculated by 



v cm ,j + v ■ 



N K 

E 

fc=i 



ran 
k 



V 

A^A 



(11) 



N h 



n 1 E( r fe- r cm,A) x K-vD 



k=l 



where II and r cm .j denote the moment-of- inertia tensor 
and the center of mass of particles in the cell, respectively. 



D. Wall Potential 

In the absence of shear flow, the monomer density pro- 
file p(y) can be calculated from the interaction potentials 
V of the dumbbells, 



p(y) = 4 



dy' e" 



p(y)/ph 



l r 

2 



erf (yK/2fc B T V 



erf ( ^K/2k B T {L y - y) 



(12) 



(13) 



where pb is the bulk monomer density, Z the partition 
function, and erf the error function. Fig. [T] shows ex- 
cellent agreement of the theoretical prediction (fl"3"|) with 
simulation data. The particles are not equally distributed 
along the wall direction; instead, at both walls, the den- 
sity is only half of the bulk density. In order to reduce 
possible slip effects, it seems desirable to make the par- 
ticle distribution as uniform as possible. An attractive 
potential is therefore applied when the center-of-mass po- 
sition of the dumbbells approaches one of the walls, 



V wa ii(ya,yi2) = -2c 2 k B T 1 



m + yi2 

2 Cl r« 



for < Cl r« ; 



Kvall(?/il,2to) 



2L y - yn - y i2 




>L y ~ Cl r^ , (14) 
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where = -^/fceT/K is the one-dimensional average 
extension of a dumbbell. The density profile is now given 



p(y) = - " dy' e~ iT^Tiy-y' f^Mv^') . ( 15 ) 

20 

The advantages of the piecewise linear form (Tl4"|) of the 

wall potential are twofold. Firstly and most importantly Q " 

it allows for an analytical integration of the equations 

of motion during the streaming step. Secondly, the den- 1 

sity profile in the absence of flow can again be calculated 

analytically (see Appendix for details). 
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FIG. 2: (Color online) Monomer density profiles at vari- 
ous dimensionless shear rates 7/wo, ranging from j/uio =0.0 
to 0.447, when attractive wall potentials are applied. The 
spring constant of dumbbells and the collision time are as 
same as those in Fig. [1] The small ripples in the profile at 
large A //uuo are due to inhomogeneities in the temperature pro- 
file, since velocity rescaling is not sufficiently efficient at high 
shear rates. The ripples to not appear for MPC-AT. 
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FIG. 1: (Color online) Monomer density profiles with 
(squares) and without (circles) attractive wall potentials ap- 
plied along the wall direction when particles approach close 
to walls. The dashed and dotted lines are the theoretical 
prediction described in Eq. (|12|l and Eq. (|15|l . respectively. 
The spring constant of dumbbells and the collision time are 
K = 0.2 and h = 0.02, respectively. Both simulations are with 
absence of shear flow. 

The simulated density profile shows excellent agree- 
ment with the analytical solution of Eqs. ([T4"|) and (fT5)) 
(see Appendix). The factors c\ and c 2 are chosen to 
obtain a nearly uniform density distribution. This is 
achieved for c\ — 1.3 and C2 = 0.4. As shown in Fig. [TJ 
the densities of point particles at both wall boundaries 
deviate by less than 10% lower from the bulk value, when 
the attractive wall potentials are applied. Simulations 
are also performed on systems of dumbbells with vari- 
ous spring constants, ranging from Ka^/feeT = 0.1 to 
Ka§/fceT = 5.0, in the absence of shear flow. It is found 
that for the given values of c\ and C2, the density pro- 
files are essentially independent of the spring constant 
of the dumbbell in the range 0.1 < Kal/k^T < 1.0. In 
Fig. [2 we plot density profiles in shear flow. At lower 
shear rates, i.e. j/ujQ < 0.1, nearly identical profiles are 
obtained as without flow. For higher shear rates, devi- 
ations of the density profile from the equilibrium profile 



become significant. Nevertheless, these profiles are still 
more uniform than those without an attractive wall po- 
tential. Our investigations are mainly focusing on rela- 
tively low shear rates, where the non-uniformity of the 
density profile is not significant. 



E. Stress Tensor and Shear Viscosity 

In the MPC model, the visc osity r? consists of a ki- 
netic and collisional contribution[33,[3y|. At steady shear 
rates, with flow along the x direction and gradient along 
the y direction, r\ is calculated by measuring the xy com- 
ponent of the stress tensor, a xy = cr^ n + er™ 1 , so that 

In the streaming step, cr^™ is proportional to the flux 
of the x momentum crossing a plane normal to the y di- 
rection. Since the stress tensor is independent of the po- 
sition of the plane, we choose y = or y = L y to measure 
the momentum transfer. In two-dimensional simulations, 

Ni 
i— 1 

where t w £ [t, t+h] is the time at which particle i bounces 
back from the wall, v x ^(t w ) and v' x ^{tw) are the velocities 
just before and after the collision with the wall, and N\ 
denotes the number of particles which hit one of the walls 
in the time interval [t, t+h]. In the collision step, particles 
close to the wall will change their velocities due to the 
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multi-particle collisions with virtual wall particles with 
average velocity v x — zL^L y , so that 



col 
>xy 



No 
i—1 



i(t + h)-v x ,i(t + h)] 



(17) 



Here N 2 denotes the number of particles which have 
multi-particle collisions with virtual particles, while 
v X} i(t + h) and v' xi (t + h) are the velocities of particle 
i before and after the collision step, respectively. In our 
simulations, N 2 is found to be much larger than N\ for 
small collision times h, indicating that the collisional part 
dominates the shear viscosity. Simulations are first per- 
formed on a system of pure point-like fluid particles to 
verify the measurement of the zero-shear viscosity from 
Eqs. (fT6|) and (fl7|) . We get p erfect agreement with the 
theoretical predictions |35l l3a \3l\ for r\. 

The shear viscosity can also be measured from system 
under Poiseuille flow[6ll38j by 



n 



pgL 2 y 

8^max 



where g is the gravitation field, and v, 
flow velocity. 



(18) 

, n;l(X the maximum 



F. Storage and Loss Moduli 

In an oscillatory shear flow, the shear rate j(t) is time- 
dependent, 



j(t) = 7o w cos (wt) , 



(19) 



where 70 and u> are the strain amplitude and the oscil- 
lation frequency, respectively. Note that the frequency 
u in Eq. (|19p is independent of the angular frequency 
luq of harmonic dumbbells in Section [TTJ In our simula- 
tions, we choose 70 1 in order to investigate the linear 
viscoelastic regime. The stress tensor is divided into two 
contributions, the viscous part a' and the elastic part a", 
so that 

<J X y(t) = a' sin (ujt) + a" cos (ut) 

= 70 [G'{uj) sin (cut) + G"(lu) cos (ujt)] , (20) 

where G' is the storage modulus, which measures the in- 
phase storage of the elastic energy, and G" is the loss 
modulus, which measures the out-of-phase energy dissi- 
pation. For a simple Maxwell fluid, G' and G" are given 
by.39] 



G' = G* 



G" = G' 



1 + (co/u>*) 2 

lu/lu* 
1 + (u/u*) 2 



(21) 
(22) 



where uj* is a characteristic relaxation frequency, and G* 
is a characteristic shear modulus. In the limit of uj <C to*, 
the loss modulus is G" — i] w, where 77 is the zero-shear 
viscosity. 



G. Kinetic Theory of Dumbbells in Solution 

In order to estimate the rheological properties of our 
model fluid, we modify the kinetic theory for dilute solu- 
tions of elastic dumbbells ■ For Hookean dumbbells 
in a solvent, the viscosity r/o, the storage modulus G' and 
the loss modulus G ' are given bv[4Cl| 



P k B T 



G'q — 



pk B T (w/w s ) 2 



G'q = rj s u + 



1 + {u/uj s f 



pk B T uo/lo s 



1 + [Oj/Ojsf ' 



where 



4K 



(23) 



(24) 



(25) 



(26) 



with solvent viscosity rj s and friction coefficient £ s of a 
monomer. Moreover, the expectation value for the square 
of the monomer separation, divided by its equilibrium 
value, is given by[ 



(r 2 ) 
(r 2 ) 



= 1 



(27) 



cq 



In Ref. |40j, the friction coefficient is obtained from 
Stokes' law for a bead of radius r in the solvent, i.e. Cs — 
677 rj s r. However, in the MPC dumbbell fluid, there ex- 
ists no explicit solvent and the monomers are point par- 
ticles instead of spheres. Nevertheless, the motion of 
the monomers is governed by the friction caused by the 
surrounding monomers which can be considered to take 
the role of the solvent. Using £ = k B T/D, which fol- 
lows from the Stokes-Einstein relation, we can thus re- 
late the friction to the diffusion constant D of a MPC 
fluid of point particles with the same monomer density. 
Similarly, we substitute the viscosity of the solvent, rj s , 
by the corresponding viscosity t^mpc of a MPC fluid of 
point particles. Theoretical expression for t/mpc an d 
D for the different collision methods can be found in 
Refs. H, HE H HO and Refs. H M, El respectively. 
The zero-shear viscosity then reads 



V = Vmpc + 



P k B T 
2 io B 



where we have introduced 

4K 4DK 



k B T 



(28) 



(29) 



Note that the limit K — > 00 corresponds to a MPC fluid 
of Np/2 point particles of mass m c . Here, the second 



term in Eq. 



vanishes, and since ^mpc(/°/2, 2m) 
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??MPc(Pi m ) for not too small p and sufficiently small h 
(so that the collisional part of the viscosity dominates), 
the viscosity resulting from this simple theory approaches 
the correct value in this limit. 

Consequently, we use the same substitutions for the 
storage and loss modulus, and for the average dumbbell 
extension, and obtain 



G' 



pk B T (uj/lo h ) 2 



and 



G" = ?7mpc^ 



(r 2 



1 + {uj/whY 



pk B T lo/ujh 



1 + (lu/ujh) 2 



(30) 



(31) 



(32) 



We emphasize that the above expressions only serve as a 
semi-quantitative description of the MPC dumbbell fluid. 
For example, the employed expressions for the diffusion 
constant neglect hydrodynamic interactions, which be- 
come important for small time steps ft. 



In Tab. HI the theoretical values for the diffusion con- 
stant D are given for ft = 0.1 for the different collision 
methods and various monomer densities [37], l4ll ]. The 
corresponding results for other time steps ft can be ob- 
tained by employing the linear relationship between D 
and ft. 



9 


£)(SHD) 


D (AT-a) 


D (AT+a) 


10 


0.1222 


0.1222 


0.1353 


20 


0.1105 


0.1105 


0.1162 


40 


0.1051 


0.1051 


0.1078 



TABLE I: Diffusion constants D of point-particle fluids for 
the standard MPC-SRD algorithm, as well as for MPC-AT-a 
and MPC-AT+a simulations for various monomer densities, 
in two dimensions. All data are calculated for collision time 
h — 0.1. Diffusion constants for other time steps h can be 
obtained by employing the linear relationship between D and 
h. Note that the values for MPC-AT— a are identical with 
those for MPC-SRD with collision angle a = 90°. 



B. Steady Shear Flow 



III. RESULTS 
A. Dimensionless Variables and Parameters 

In the remainder of this article, we introduce dimen- 
sionless quantities by measuring length in unit of the 
lattice constant do, mass in unit of the dumbbell mass 
m c , time in units of ao\/m c /k^T, velocity in units of 
\J fceT/m , monomer number density p in units of d , 
where d is the spatial dimension, and the spring constant 
K in units of k-QT/a^. The shear rate 7 and all kinds 
of frequencies are measured in units of \/k^T /m c ag. Fi- 
nally the viscosity r\ is in units of y 'iii^bT / 'a^, and the 
storage modulus G' and the loss modulus G" are in units 
of k^T/a^. In these dimensionless units, the mean free 
path A (in units of the lattice constant) becomes equiva- 
lent to the time step ft. 

In our simulations, harmonic dumbbells with N p 
point particles are initially placed in a two- or three- 
dimensional rectangular box at random. We choose the 
average number density of point particles p — 20 and 
L x = 50 for all two-dimensional simulations which re- 
sults in N p = lOOOLj,. The collision time ranges from 
ft = 0.01 to ft = 0.2, while the spring constant ranges 
from K = 0.1 to K = 5.0. The rotational angle is cho- 
sen a = 90° and a = 130° for two- or three-dimensional 
simulations, respectively. We use small ft and large a 
to obtain large Schmidt numbers required for fluid-like 
behavior [37], EH- Most of the results shown are ob- 
tained from two-dimensional systems, except in a few 
cases where this is explicitly mentioned. 




FIG. 3: (Color online) Snapshots of dumbbell configurations 
in steady shear flow. The system size is L x — L y = 50. 
Half of each dumbbell is colored red, the other half yellow for 
reason of visualization. In each frame only 2500 dumbbells 
are shown, so that the density is 10 times as high as appears 
from the pictures. The spring constant and the collision time 
are K = 0.2 and h = 0.02, respectively. From (A) to (D), the 
applied shear rates are 'j/ujh ~ 0.0565, 0.565, 2.83 and 5.65. 

In Fig. [31 we present snapshots for steady shear flow 
with a simulation box containing 25000 dumbbells. At 
lower shear rates, i.e. a //ujh < 0.6, see Fig. [3^ and[3f3, 
the average extension of the dumbbells is hardly dis- 
tinguishable from the equilibrium value. In these two 
cases, the shear flow is not strong enough to align the 
dumbbells along the flow direction, so that both systems 
are still isotropic. With increasing 7, shear forces over- 
whelm entropic forces. As a result, dumbbells are largely 
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stretched, at the same time reorientated along the flow 
direction, as presented in Fig.[3p and[3p. Note that near 
both the walls, the average size (r 2 ) 1 / 2 of the dumbbells 
in flow is larger than in the bulk. Also, an alignment 
of the dumbbells is found near the walls, both with and 
without shear flow, with peaks at y = and y = L y . 
This is an effect of the geometrical constraints imposed 
on anisotropic particles by a hard wall. Furthermore, 
a maximum of the extension occurs at a finite distance 
from the wall, which we attribute to the combined ef- 
fect of the wall and the flow conditions; dumbbells very 
close to the wall are sterically oriented parallel to the wall 
and thus experience only a very small shear force, while 
those a little further away are close to the average in- 
clination angle (see Fig. IIII Bl below), which corresponds 
to the largest stretching. The distance of the position 
of the maximum from the wall decreases with increasing 
shear rate, and seems to approach the size of the colli- 
sion cells for large 7. The relative peak height increases 
with increasing shear rate. For example, we find that the 
maximum extension (r 2 ) 1 / 2 near the wall is about 11% 
larger than the bulk extension for 7 /ujh = 1.13, while it 
is about 28% larger than in the bulk for j/u>H — 2.83. 



angle is defined as the angle between the average orien- 
tation of the end-to-end vector of a dumbbell and the 
the flow direction. At lower shear rates, j/u>H < F the 
inclination angle approaches 9 = 45°, while it decays to 
zero for large shear rates with a power law 7 . 
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FIG. 4: (Color online) Distribution of dumbbell configura- 
tions for the system shown in Fig. [3] Each dot indicates the 
end-to-end vector of a dumbbell. 



FIG. 5: (Color online) The inclination angle 9 as a function 
of dimensionless shear rate ^/uh for the system of Fig. [3] 
with spring constant K = 0.2, collision time h = 0.02, and 
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Fig. IIIIBI presents the extensional and orientational 
distribution of the dumbbells for various shear rates. At 
lower shear rates, 7/wy < 1, the end-to-end vector of 
the dumbbells is distributed on a circle, see Fig. IIII BA 
and B, indicating an isotropic orientation. At a higher 
shear rate, j/uH — 2.83, the orientational distribution 
becomes an elongated ellipse, see Fig. IIIIBC . With in- 
creasing shear rate, the distribution elongates further. 
Simultaneously, dumbbells become more aligned with the 
flow direction, as can be seen quantitatively from the in- 
clination angle 9 shown in Fig. [5J Here, the inclination 



FIG. 6: (Color online) Shear viscosity rj and scaled average 
dumbbell length (r 2 ) 1//2 /r{, 2 ' as a function of dimensionless 
shear rate j/uiH- Systems with the wall separation L y = 
10, 20, 30, and 50 are investigated. The spring constant is 
K = 0.2 and the collision time h — 0.02. 

In Fig. [6l we plot the shear viscosity 77 as a function 
of dimensionless shear rate j/uJH for various wall sep- 
arations L y ranging from 10 to 50. In each system, r] 
remains constant until the applied shear rate reaches a 
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critical value jc/uiH ~ 5. The shear viscosity then decays 
rapidly as 7 further increases, showing a typical "shear- 
thinning" behavior. Fig. [6] also shows the average exten- 
sion of dumbbells (r 2 ) 1 / 2 /r[ ) 2 ' as a function of the shear 
rate. 

Two comments are required here. Firstly, in our MPC 
model, an entanglement between dumbbells is not taken 
into account, so that they can freely cross each other. 
Also, the absence of an excluded-volume interaction im- 
plies that there is no benefit of a para-nematic order- 
ing in terms of an increased sliding of parallel dumb- 
bells along each other as in solutions of rod-like col- 
loids; instead, parallel dumbbells interact very similarly 
to isotropically oriented dumbbells, since in both cases 
the monomers colloide with other monomers in exactly 
the same fashion. Thus, our system is very similar to a so- 
lution of non- interacting harmonic dumbbells, for which 
- in the absence of a finite extensibility - neither "shear- 
thinning" nor "shear-thickening" is expected[3l.l40j. com- 
pare Eq. (j2li|) . Secondly, the size of the simulation box 
should have no influence on the bulk viscosity at a given 
shear rate. However, the plateau value of the viscosity 
increases strongly with the wall separation L y . This in- 
dicates that boundary effects could be responsible for the 
observed "shear-thinning" behavior. 
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little slip at the walls, as expected. However, the velocity 
decays rapidly in a boundary layer of thickness A, then 
decays linearly to zero at the middle plane. Obviously 
the applied shear rate 7 is not appropriate to calculate 
the shear viscosities from the stress tensor a xy by 77 = 
a xy /j. An effective shear rate j e g is therefore introduced 
instead, which characterizes the linear bulk part of the 
velocity profile. At a given shear rate, the larger the wall 
separation, the less the effective shear rate deviates from 
7, since the finite-size effect is much stronger in smaller 
systems. The ratio j/%s between the applied and the 
effective shear rates is plotted in Fig. [HJ as a function 
of j/u>H- At lower shear rates, i.e. 7 < j c , where % 
is the critical shear rate, the ratio j/jes is independent 
of the shear rate. When the applied shear rate becomes 
larger than this critical value, the effective shear rate 7 e ff 
increases more slowly than 7. 
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FIG. 8: (Color online) Ratios between the applied shear rates 
7 and the effective shear rates j e g as a function of a //ljh for 
the same systems as in Fig. [6] 

Consequently, the effective shear viscosity can be cal- 
culated by 



7cff 



(33) 



FIG. 7: (Color online) Velocity profiles for wall separations 
L y = 10, 20, 30, and 50. Data are obtained for spring 
constant K = 0.2, collision time h — 0.02, and shear rate 
"//ujh = 0.565. Fhe dashed line corresponds to the applied 
shear rate. Solid lines represent fits to the bulk part of the 
velocity profiles, their slopes yield the effective shear rates 
7cff- 

We therefore examine the velocity profiles for systems 
with various wall separations L y . In Fig. [7J the average 
velocities v x of the monomers along the flow direction 
are plotted as function of L y for a fixed shear rate of 
j/wji = 0.565. The velocities at the boundaries deviate 
only very little from the wall velocities, i.e. there is very 



In Fig. [51 rjes is shown against %g /u>h for various wall 
separations L y . The data for different system sizes now 
all fall onto a single master curve, which describes the 
bulk shear viscosity. Now, instead of "shear thinning" 
shown in Fig. a very weak "shear thickening" behavior 
is observed when 7 e ff/wH > 1. Three-dimensional simu- 
lations are also carried out for systems of 20 x 20 x 10 
boxes along the x, y, and z directions. For the same pa- 
rameters K = 0.2 and h — 0.02, weak "shear thickening" 
behavior is also observed, as shown in the inset of Fig. [9l 
when jeg/ujH reaches the critical value, 7 c ,cff/ w ff ~ 2. 
The value of the critical shear rates are found to be very 
similar in two and three dimensions. 
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FIG. 9: (Color online) Master curve of the viscosity rj e f; as a 
function of the effective shear rate 7 e g on a semi-logarithmic 
scale. Symbols are the same as in Fig. [6] and Fig. [8] For 
strong shear flow, i.e. j e ff/u)H > 2, the viscosity increases 
("shear thickening"). The dashed line is fitted to the data 
with wall separations L y > 20. In the inset, r/ e ff is plotted 
as a function of A j c s/^H for three-dimensional systems. A 
20 x 20 x 10 simulation box is chosen, while the spring constant 
K, the collision time h and the average number density p are 
the same as in the two-dimensional systems. 



Fig. [5] shows that the effective shear viscosity r/ c g 
is nearly independent of the shear rate for jog/LUn < 
jc,es/i^H ~ 2. This critical shear rate corresponds to the 
onset of the apparent "shear thinning" observed in Fig. [51 
as well as the deviation of j/jes from its low-shear-rate 
value in Fig. [5] It should be noticed that the value of 
iccs/^H ~ 2 implies jc/^H is in the range [3,6.4] for 
system sizes L y € [10,50], compare Fig. [5) However, it 
is important to note that there is already a pronounced 
alignment and stretching of the dumbbells for smaller 
shear rates; Fig. [5] shows that the inclination angle 9 
has decreased from 9 = 45° in the absence of shear flow 
to 9 sa 15° at j/ujh — 3, while Fig. [5] indicates that 
( r 2)i/2/ r (2) w 2 at A//uJh = 3 

The spring constant K of the dumbbells is of great 

importance, since it controls the elasticity of the fluid. 

We have therefore examined velocity profiles of systems 

of dumbbells with various spring constants. In Fig. [TUl 

the simulation results are plotted for a fixed applied shear 

rate 7 = 0.01. The effect of the boundary layer becomes 

more pronounced with decreasing spring constant. By 

fitting the linear parts of the velocity profiles, we find 

that, for the same shear rate 7 = 0.01, the effective shear 

rate j c g for dumbbells with K = 0.1 is about 10 times 

lower than that with the highest spring constant studied 

here, K = 4.0. The thickness of the boundary layer is 

. (2) 

proportional to the equilibrium average extension Tq = 
yj2k^TfK, as shown in the inset of Fig. [TUl 



FIG. 10: (Color online) Velocity profiles for various spring 
constants, ranging from K = 0.1 to K = 4.0. The wall sepa- 
ration in each case is L y = 10. The dashed line corresponds 
to the applied shear rate 7 = 0.01, while solid lines are the 
fitted effective velocity profiles. The inset shows the thickness 
of the boundary layer A as a function of equilibrium extension 
rg 2 ' = ^/2fcBT/K. The dashed line is a linear fit. 



The zero-shear viscosity ?y c ff is found to depend linearly 
on 1/K, see Fig.QTj As K increases, the effective viscos- 
ity VeS approaches the expected value of system of point 
particles with mass of m c and density p/2. The same lin- 
ear relationship between r] c s and 1 /K is also obtained in 
three-dimensional systems, as shown in Fig. 1111 Not only 
the linear dependence of r) e g on 1/K but also the pref- 
actors are in very good agreement with the theoretical 
predictions |J5SJ. 

The scaled mean free path A, which determines how 
far a point particle travels between collisions, is another 
important parameter which affects the shear viscosity. 

We always employ small mean free paths mm, so that 

the collisional viscosity is dominant compared to the ki- 
netic viscosity. The data of Fig. [T2Ta) demonstrate that 
the zero-shear viscosity increases linearly with 1/A, for all 
spring constants K studied here, as it does for a system of 
point particles [H, [H, H3]- However, the slope decreases 
with increasing K, in good agreement with the analytical 
results obtained from Eq. |28|) . as shown in the inset of 
Fig. [Ha). 

The weak "shear-thickening" behavior is observed for 
all mean free paths investigated here, see Fig. [T2] Thus, 
this weak "shear-thickening" behavior is intrinsic to the 
MPC algorithm, and cannot be avoided by a variation 
of the collision time. Fig. [TWa) indicates that the crit- 
ical shear rate 7 c ,eff/wo depends only very weakly on 
the mean free path A. Therefore, we present the sim- 
ulation data in Fig. [T^] as a function of A / e f{/ujQ 1 since 
luq = (2K/m) 1 / 2 is independent of A, while ujh decreases 
linearly with A. The "shear-thickening" behavior be- 
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FIG. 11: (Color online) The zero-shear viscosity r) c s as a 
function of spring constant K in both two-dimensional (cir- 
cles) and three-dimensional (squares) systems. The solid 
lines are linear fits, the dashed lines indicate the theoreti- 
cal predictions (|28[1 . In all simulations, the collision time is 
h — 0.02. Two-dimensional simulations are performed in sys- 
tems of 50 x 10 boxes, while three-dimensional simulations 
are in 30 x 30 x 20 boxes along the x, y and z directions, re- 
spectively. The average number density in three-dimensional 
systems is p = 10, which is half of value in two-dimensional 
systems. 



comes more pronounced and slowly shifts to smaller val- 
ues of 7 e ff /u>o for system of dumbbells with smaller spring 
constants, see Fig. H2T b). In the range of investigated 
spring constants and mean free paths, the shear thick- 
ening occurs roughly at 7 c , c ff/^o — 0.1. It is important 
to note that the viscosity of the standard point-particle 
MPC fluid is also not independent of the shear rate, but 
shows a weak shear-thinning behavior at high shear rates 
[35l | . For our model parameters and in two dimensions, 
this shear-thinning behavior sets in at a shear rate j c ~ 1. 
Thus, with increasing K, shear thickening occurs at a 
slowly increasing 7 Cje ff/wo for K < 1; for larger spring 
constants K > 5, shear thinning is observed instead, and 
7 C) eff/wo decreases again (since 7 Cje ff — * 1 and c^o — » 00 
for K — > 00). 
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FIG. 12: (Color online) The effective shear viscosity 7j e g as 
a function of the dimensionless effective shear rate 7 e ff/o;o, 
on a double-logarithmic scale, (a) For fixed spring constant 
K = 0.2 and various mean free paths A = 0.01, 0.02, 0.05, 0.1 
and 0.2. (b) For fixed mean free path A = 0.02 and various 
spring constants K = 0.2, 0.5, 1.0 and 5.0. In both cases, 
the wall separation is L y — 10. The inset in (a) shows the 
zero-shear viscosity rfeff as a function of 1/A. The solid line 
indicates the theoretical result for K — > 00, while the other 
lines show the predictions (|28[) for K = 0.2, 0.5, and 5.0. 



C. Small-amplitude Oscillatory Shear Flow 

Another way to explore the viscoelastic properties of a 
fluid is to apply a small-amplitude oscillatory shear flow. 
We use here the strain amplitudes 70 =7/0; in the range 
0.1 to 0.5 to mimic a small amplitude shearing. The fre- 
quencies u> ranges from 10~ 4 to 10" 1 in our simulations, 
which provides a wide range of shear rate from 10~ 5 to 
5 x 10~ 2 . 

The storage and loss moduli as a function of oscilla- 
tion frequency are plotted in Fig. 1131 Similarly to the 



simulations of steady shear flow, effective shear rates are 
measured from the bulk velocity profiles at times when 
cos(ujt) = ±1. By doing so, all the simulation data 
fall onto master curves at various wall separations from 
L y — 10 to 50. As can been seen from Fig. [T37 a). the 
storage modulus G' is well fitted by Eq. ((2Tj) . indicat- 
ing that the dumbbell system exhibits a typical behav- 
ior of a Maxwell fluid. The relaxation frequency uj* ob- 
tained from the fit of the storage modulus G' against u> 
is then used in Eq. |22|) to fit the loss modulus G" . In 
Fig. [T3T b). at low frequencies, u> < 0.02, the simulation 
data follow the expected linear w-dependence very well. 
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In this linear regime, the shear viscosity is then calcu- 
lated by r\ = G"(lu)/ui, which yields n = 565, in excellent 
agreement with the result in steady shear flow, compare 
Fig. [HI Note that the fitted values for the amplitude G* 
in Eqs. (JSTJ and ([22]) differ by about a factor 2. This in- 
dicates that the system investigated here does not behave 
exactly like a simple Maxwell fluid. 
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FIG. 13: (Color online) (a) Storage G' and (b) loss modu- 
lus G" , as a function of oscillation frequency w on a double- 
logarithmic scale, for systems with various wall separations 
ranging from L y = 10 to 50. The spring constant and the 
collision time are K = 0.2 and h = 0.02, respectively. The 
dashed line in (a) is fitted by the Maxwell model, Eq. (I21|l . 
on the basis of all simulation data, while the one shown in (b) 
is based on data for oscillation frequencies uj < 0.02. 

In Fig. [TU we examine the storage and loss moduli of 
system of dumbbells with various spring constants. As 
in Fig. [13j simulation results are all well fitted by the 
Maxwell equations ([21]) and ([22]) . except for somewhat 
different amplitudes G* . The relaxation frequency to* is 
found to agree very well with uih, as shown in the inset 
of Fig. [Ml At lower oscillation frequency in Fig. RUT b). 
the viscosities calculated from G"(uj)/w are 77 = 565, 253 
and 144 for systems with K = 0.2, 0.5 and 1.0, respec- 
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FIG. 14: (Color online) (a) The storage G' and (b) the 
loss moduli G" , as function of oscillation frequency u on a 
double-logarithmic scale, for systems of dumbbells with var- 
ious spring constants ranging from K = 0.2 to K = 1.0. 
The wall separation and the collision time are L y — 10 and 
h = 0.02, respectively. The inset shows the fitted relaxation 
frequencies to* as a function of the frequency ujh predicted by 
Eq. (|29p . The dashed line shows the identity u* — ujh- 



tively. These values are again in excellent agreement with 
those calculated from Eq. (|33[) in steady shear flow. For 
all spring constants K, the fitted amplitudes G* for the 
storage moduli G' are about half of those calculated for 
the loss moduli G" . This indicates that even for a system 
of dumbbell with high spring constant, a simple Maxwell 
model is not appropriate for a quantitative description. 



D. Angular Momentum Conservation 

The viscosity of a simple MPC-AT+a fluid (with 
angular-momentum conservation) is about a factor 1/2 
smaller than of a MPC-AT-a fluid HIH. We thus ex- 
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pect the viscosity of the dumbbell fluid to be affected by 
angular-momentum conservation as well. The simulation 
results for both MPC-AT-a and MPC-AT+a methods 
are compared in Fig. [T5] We find that the effective zero- 
shear viscosity t] c q increases linearly with the monomer 
density p for p > 5. The corresponding theoretical results 
(|28|) are in good agreement with the simulation results 
for both investigated spring constants. Minor deviations 
from the linear relationship of r] c g with p originate from 
the variation of the diffusion constant at low densities, 
which approaches a constant value for high p. The vis- 
cosity of MPC-AT+a is lower than for MPC-AT-a, al- 
though this effect is less pronounced than for pure point- 
particle MPC fluids, since the main contribution to the 
viscosity originates from the spring tension. 
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FIG. 15: (Color online) Effective shear viscosities Tfeg as a 
function of the density p for MPC-AT-a and MPC-AT+a, 
each for spring constants K = 0.2 and K = 0.4. The lines 
represent the theoretical results obtained from Eq. (12811 . The 
wall separation and the collision time are L y — 20 and h = 
0.014, respectively. 

In Fig. 1161 we present the average squared dumbbell 
extension, determined in the bulk as a function of the 
effective shear rate, along with the theoretical results 
(|52"j) . for both the MPC-AT-a and MPC-AT+a methods. 
Note that the diffusion constant D in Eq. (|32|) is different 
for angular-momentum conserving and non-conserving 
methods. The angular-momentum conservation has only 
little effect on the spring extension; for a given effective 
shear rate, the extension is slightly lower for the angular- 
momentum conserving method. The agreement of the 
simulation data with the theoretical result (|3"2")l is again 
remarkably good. 

IV. DISCUSSION 

As a further test for the correct calculation of the ef- 
fective viscosity by the procedure described in Sees. Ill El 



FIG. 16: (Color online) Scaled average of the dumbbell ex- 
tension, (r 2 ) / (r' 2 ) cq — 1, as a function of the effective shear 
rate ^ch/loh for angular-momentum conserving and non- 
conserving methods. The spring constant and collision time 
are K = 0.2 and time step h — 0.014, respectively, the density 
is p — 10 and the wall separation is L y = 20. The dashed line 
represents the theoretical result (I32[) . 



and IIII Bl we have also determined the viscosity from 
Poiseuille flow. As in Ref. l33l . we apply a gravita- 
tional force of strength g parallel to the walls, with g 
in the range from g = 0.0001 to g = 0.01 (in units of 
/cbT/oo). We fit the central part of the velocity profile 
to a parabolic flow curve. The value of this curve at the 
wall positions determines the effective wall slip. When 
this slip velocity is subtracted, Eq. JT5|> in Sec. Ill El is 
employed to determine the viscosity [4J|. We have used 
this method for a system of dumbbells with K = 0.2 in a 
30 x 30 box. Excellent agreement between the two meth- 
ods to calculate the zero-shear viscosity is obtained. 

Our results for the dependence of the inclination an- 
gle 9 on the shear rate can be compared with the decay 
of the inclination angle of flexible and semi-flexible poly- 
mers. For dilute polymer solutions in the asymptotic 
regime of high shear rates (where the finite extensibility 
is important), 6 has been predicted from Brownian dy- 
namics simulations (45| and theory [46| to decay with a 
power law j~ - 3 and -"y^ 1 / 3 , respectively. For extensible 
dumbbells, the theory of Ref. predicts [4^ 8 ~ 7 _1 , in 
excellent agreement with our simulation results. 

The wall slip in polymer melts has been studied exten- 
sively. In this case, molecular dynamics simulations of 
polymer fluids with Lennard- Jones interactions between 
monomers give a wall slip with a boundary layer thick- 
ness, which is on the order of the monomer diameter a 
or less 0, HU . Our model could be compared more eas- 
ily with results for polymer solutions, because our model 
does not include excluded- volume interactions. However, 
there is little knowledge about semi-dilute polymer solu- 
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tions near a wall under flow conditions. Nevertheless, 
some comparisons with polymer melts with moderate 
chain lengths are possible, where entanglement effects 
are absent. For example, the molecular dynamics sim- 
ulations of Zhang et al. [44j| show a maximum of mean 
squared radius of gyration at a finite distance A m from 
the wall, which shifts from A m ~ 1.5a for chains with 4 
monomers to A m ~ 2.2a for 10 monomers. 



V. SUMMARY AND CONCLUSIONS 

A multi-particle collision dynamics (MPC) algorithm 
has been developed to investigate the viscoelastic prop- 
erties of harmonic-dumbbell fluid in shear flow. The 
method is based on alternating streaming and collision 
steps, just as the original MPC method for Newtonian 
fluids. The only modification is to replace the ballistic 
motion of fluid point particles by harmonic oscillations 
during the streaming step. In this model, the entan- 
glement between dumbbells is neglected. Moreover, the 
storage and loss moduli are calculated by introducing a 
small amplitude oscillatory shear flow. 

Our results can be summarized as follows: 

First, under steady shear flow, the dumbbells keep 
their isotropic distribution at low shear rates, but get 
highly stretched and orientated along the flow direction 
at high shear rates. The velocity profile is not uniform 
along the gradient direction, but boundary layers of high 
shear develop near the walls. The thickness of these 
boundary layers is found to scale with the size of the 
dumbbells in the absence of flow. The effective shear vis- 
cosity, calculated from the ratio between the off-diagonal 
component of the stress tenor, a xyi and the effective shear 
rate 7 e ff , expresses a very weak "shear thickening" behav- 
ior at high shear rates. 

Second, the dependence of the viscosity on two pa- 
rameters, the spring constant K of the dumbbells and 
the collision time h, has been investigated. These two 
parameters are of central importance, since the former 
controls the elastic energy of the system, while the lat- 
ter determines the mean free path A, which measures the 
fraction a the cell size that a fluid particle travels on aver- 
age between collisions. We find that the shear viscosity of 
the dumbbell fluid increases linearly with 1/K and with 
h. 

Third, the storage and loss moduli of our viscoelastic 



solvent are studied by imposing an oscillatory velocity on 
the two solid walls. The storage modulus G' is found to 
be proportional to u 2 at low frequencies, and to level off 
at to* — u>h- Its behavior over the whole frequency range 
studied here is well described by a Maxwell fluid. The loss 
modulus G" increases linearly with lu for low frequencies. 
The shear viscosities obtained from the ratio G" /u) at 
low shear rates agree very well with those obtained from 
simulations with steady shear. On the other hand, for 
lu > loh, we find that the data approach a plateau value, 
while for a Maxwell fluid G" would decrease again for 
higher frequencies. 

Our numerical results are quantitatively in good agree- 
ment with a simple theory, based on the kinetic theory 
of dilute solutions of dumbbells, where the transport co- 
efficients of the standard MPC point-particle fluid are 
employed for the viscosity and the diffusion constant of 
the solvent. 

In our MPC algorithm of harmonic dumbbells, both 
elastic and viscous behaviors of solvent particles can be 
modeled properly, while hydrodynamic interactions are 
efficiently taken into account. These are valuable assets 
to guide future simulations on investigating rheological 
properties of suspensions of spherical, rod-like or poly- 
meric solute molecules in viscoelastic fluids. 
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Appendix: Analytical solution of the density profile 
with attractive wall potentials 

Combining Equations (| 14[) and (|15[) . the density pro- 
file, when attractive wall potentials are introduced, can 
be solved analytically. Considering the symmetry of the 
density profile, p(y) — p(L y — y), only the initial half part 

need to be taken into account. For < y < 2c\v^\ we 
then arrive at 



p{y) = - { / dy> e-^v-v'f/^T 



2c ir ( 1} -y 



dy e 



' -K(y-y') 2 /2k B T 



5 2c 2 k B T(l-(y+y')/2c 1 r ( 1) ) _ 1 



14 



which implies 

P(V) = \ {erf ( V / K/2fc B T (L y - y)) + erf (^K/2k B T (2y - 2c 1 r< 1) )) 

-c 2 k B T + cir^yK/fcer i 
Cl r^V2K/fc B r / 

I 

while for 2c]Tq < y < L y /2, the density profile is given 
byEq. (H3D. 



exp 



5 Tr«) 2 



4c 2 fc B T(c4 1) - 



2 Cl r« 
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